Data Prep and Display
source("ozone_krige.R")
# get elevation from DEM data.
# 1. Ozone with Temperature
# 2. Ozone with distance to nearest A1 road
# 3. Ozone with total length of all A1 roads within 500 m
# 4. Ozone with temperature and RH
# 5. Ozone with temperature and distance to nearest A1 road
# 6. Ozone with temperature and elevation
# 7. Ozone with temperature and elevation and distance to nearest A1 road.
path_to_data_folder="../data/krige_data/"
added_tiffs_list = list.files(path_to_data_folder)
# Adding Max Relative Humidity
o3_data_to_add = raster(paste0(path_to_data_folder,added_tiffs_list[8]))
added_data_projected = raster::projectRaster(o3_data_to_add, crs=prg)
o3_projected$rh=round(raster::extract(added_data_projected,o3_projected),2)
# Adding Temperature
o3_data_to_add = raster(paste0(path_to_data_folder,added_tiffs_list[14]))
added_data_projected1 = raster::projectRaster(o3_data_to_add, crs=prg)
o3_projected$max_temp=round(raster::extract(added_data_projected1,o3_projected),2)
plot(added_data_projected, main= "Maximum Relative Humidity")
plot(den_projected, col="red", add=TRUE)
plot(o3_projected, col="green", add=TRUE)
plot(o3_500m, col="blue", add=TRUE)
plot(added_data_projected1, main="Maximum Temperature")
plot(den_projected, col="red", add=TRUE)
plot(o3_projected, col="green", add=TRUE)
plot(o3_500m, col="blue", add=TRUE)
head(cbind(o3_projected$site_name,o3_projected$rh,o3_projected$elev))
## [,1] [,2]
## [1,] "Aurora East" "71.64"
## [2,] "Boulder Reservoir" "62.7"
## [3,] "Denver - Camp" "69.21"
## [4,] "Highland Reservoir" "70.22"
## [5,] "La Casa" "69.21"
## [6,] "National Renewable Energy Labs - Nrel" "66.11"
# join elevation column back to o3_projected
# o3_projected = merge(x=o3_projected,y=o3_elevs,by="site_name", duplicateGeoms = TRUE)
#o3_projected
year_o3 = as.data.frame(o3_projected)
year_o3 = year_o3 %>%
dplyr::select(c("site_name","max_temp","rh","long","lat"),everything())
# use to create dataframe of specific months, ex below is summer
# summer_o3 = year_o3 %>%
# dplyr::select(contains(c("site_name","lat","long","elevation","Apr","May","Jun","Jul","Aug","Sep","Oct")))
# preview data
#summer_o3
Krige and Error Results
# # create empty list to store information
# kriges = vector("list")
# kriges.cv = vector("list")
# rmses = vector("list")
# preds = vector("list")
# vars = vector("list")
#
# # loop thru columns
# # reformulate to create formula based on var_name for each month
# # krige + store
# # krige.cv and rmse + store
# for(i in seq_along(cols)) {
# var_name = cols[i]
# k = autoKrige(reformulate("max_temp+rh",var_name),year_o3,grd, model = c("Ste","Sph","Mat","Exp","Nug","Gau","Lin"))
# k.cv = autoKrige.cv(reformulate("max_temp+rh",var_name),year_o3,model = c("Ste","Sph","Mat","Exp","Nug","Gau","Lin"))
# rmse = sqrt(mean(k.cv$krige.cv_output$residual^2))
# mean_mth_pred = mean(k$krige_output$var1.pred)
# mean_mth_var = mean(k$krige_output$var1.var)
# kriges[[i]] = k
# kriges.cv[[i]] = k.cv
# rmses[i] = rmse
# preds[i] = mean_mth_pred
# vars[i] = mean_mth_var
# }
#
# # plot
# for(i in seq_along(cols)) {
# month = names(year_o3[,(2+i)])
# plot(kriges[[i]], ylab = "Semivariance", xlab = month)
# }
#
# # add RMSE, predictions, and variance to gt table
# rmse_tibble = tibble(month_year = cols,
# rmse = rmses,
# mean_predicted_o3 = preds,
# mean_variance_o3 = vars)
# rmse_tibble
# rmse_gt = gt(rmse_tibble)
# rmse_gt
